EREG is the core onco-immunological biomarker of cuproptosis and mediates the cross-talk between VEGF and CD99 signaling in glioblastoma

Background Glioma is the most prevalent primary tumor of the central nervous system. Glioblastoma multiforme (GBM) is the most malignant form of glioma with an extremely poor prognosis. A novel, regulated cell death form of copper-induced cell death called “cuproptosis” provides a new prospect for cancer treatment by regulating cuproptosis. Methods Data from bulk RNA sequencing (RNA-seq) analysis (The Cancer Genome Atlas cohort and Chinese Glioma Genome Atlas cohort) and single cell RNA-seq (scRNA-seq) analysis were integrated to reveal their relationships. A scoring system was constructed according to the cuproptosis-related gene expression, and core genes were experimentally verified using real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR), Western blot (WB), immunohistochemistry (IHC), and immunofluorescence (IF). Moreover, cell counting kit-8 (CCK8), colony formation, 5-ethynyl-2’-deoxyuridine (EdU) incorporation, transwell, and flow cytometry cell cycle were performed to evaluate cell proliferation, invasion, and migration. Results The Cuproptosis Activation Scoring (CuAS) Model has stable and independent prognostic efficacy, as verified by two CGGA datasets. Epiregulin (EREG), the gene of the model has the most contributions in the principal component analysis (PCA), is an onco-immunological gene that can affect immunity by influencing the expression of programmed death-ligand 1 (PD-L1) and mediate the process of cuproptosis by influencing the expression of ferredoxin 1 (FDX1). Single cell transcriptome analysis revealed that high CuAS GBM cells are found in vascular endothelial growth factor A (VEGFA) + malignant cells. Oligodendrocyte transcription factor 1 (OLIG1) + malignant is the original clone, and VEGF and CD99 are the differential pathways of specific cell communication between the high and low CuAS groups. This was also demonstrated by immunofluorescence in the tissue sections. Furthermore, CuAS has therapeutic potential for immunotherapy, and we predict that many drugs (methotrexate, NU7441, KU -0063794, GDC-0941, cabozantinib, and NVP-BEZ235) may be used in patients with high CuAS. Conclusion EREG is the core onco-immunological biomarker of CuAS and modulates the cross-talk between VEGF and CD99 signaling in glioblastoma, and CuAS may provide support for immunotherapy and chemotherapy. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-023-03883-4.


Background
Glioma is the most prevalent primary tumor of the central nervous system [1]. It is classified into grades I to IV by the World Health Organization based on its malignant features, wherein grades I, II, and III are lowgrade glioma and grade IV is high-grade glioma, also known as glioblastoma multiforme (GBM) [1,2]. GBM is the most malignant form of glioma with an extremely poor prognosis. Despite considerable advances in the development of treatments, including surgical resection, radiotherapy, and chemotherapy, little progress toward prolonged survival and better prognosis has been achieved over the last few decades [3]. The modest median overall survival (OS) time in GBM is approximately 14 months, and only 5% to 6.8% of patients with GBM survive 5 years after diagnosis [1][2][3][4]. Multiple clinical trials, including those on immunotherapy, have been conducted for patients with GBM; however, the results did not conclude the expected results [1][2][3][4]. In our previous study, we generated a ferroptosis-related prognostic risk score model to predict the clinical significance and immunogenic characteristics of GBM [5]. However, the biomarkers and predictors for patient outcome and the immunotherapy response of GBMs have not been fully elucidated, and existing predictive models are far from satisfactory.
Beyond classical apoptosis, several forms of regulated cell death (RCD) have been identified, such as ferroptosis, necroptosis and pyroptosis [6]. These RCD subroutines differ in the initiating stimuli, intermediate activation events, and end effectors. Various heavy metals are essential micronutrients; however, the insufficiency or excessive abundance of metals can trigger cell death, which can induce RCD through different subroutines. For example, ferroptosis has been defined as an iron-dependent form of oxidative cell death caused by unrestricted lipid peroxidation [7]. A novel RCD form of copper-induced cell death called "cuproptosis" was proposed by Tsvetkov et al. [8], which is gaining attention in the field. Cuproptosis differs from oxidative stress-related cell death (e.g., apoptosis, ferroptosis, and necroptosis). In contrast, mitochondrial stress, especially the aggregation of lipoylated proteins and destabilization of Fe-S cluster proteins, results in proteotoxic stress and ultimately cell death. Hence, it may provide a new prospective for cancer treatment by regulating cuproptosis.

Cuproptosis activation scoring model
Ten cuproptosis characteristic genes (FDX1, LIAS, LIPTI, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) were obtained from previous studies [8]. Based on these characteristic genes, consistent clustering analysis was used to identify two sample clusters in the TCGA-GBM sample. The gene expression count data of two sample clusters were calculated based on R packet DESeq2 to identify differentially expressed genes (DEGs) [10]. Candidate prognostic genes were identified from differentially expressed genes based on a univariate Cox regression analysis, and redundant factors were further filtered using LASSO Cox analysis. Based on the contribution of the prognostic genes to principal components 1 and 2, CuAS was defined as: where HR (hazard ratio) is derived from the Cox analysis. The CuAS model in the validation cohort was reproduced with a similar formula. The construction of the cluster and CuAS models is illustrated in a schematic diagram (Fig. 14).

Biofunction prediction
The GO/KEGG enrichment analysis based on GSVA calculated the GO/KEGG signature activity scores in cuproptosis activation score(CuAS) = Gene HR>1 * (PC1 + PC2) − Gene HR<1 * (PC1 + PC2) The normal human astrocyte cell line HA1800 and human glioma tumor cell lines U87, U251, LN229 and A172 were purchased from the Cell Bank of the Chinese Academy of Sciences. The STR identification reports of the cancer cell lines are presented in Additional materials (see Additional materials-cell lines STR identification), and we also used CCLA, an excellent cell line identification database, for secondary identification to ensure no cross-contamination of cell lines [15]. The cells were cultured in Dulbecco's Modified Eagle's Medium (DMEM) (Gibco) containing 10% heat-inactivated fetal bovine serum and 1% penicillin/streptomycin. qRT-PCR was conducted to compare the gene expression in 20 tumor samples in adjacent normal tissue. qRT-PCR was performed in triplicate using samples derived from three independent experiments. Formalin-fixed, paraffinembedded GBM tissues were used for IHC staining. The primers' sequence are provided (see Additional file 1: Table S3).

Lentivirus infection assay
The assay complies with the protocol described in a previous article [16]. Short hairpin RNA (shRNA) against EREG (shEREG) and a negative control shRNA (sh NC) were designed and synthesized by GeneCreate (Wuhan, China). The lentivirus, pLent-shEREG-Flag-Puro or its negative control (NC) pLent-Flag-Puro (GeneCreate) was used to infect GBM cells with enhanced infection solution (GeneCreate) according to the manufacturer's protocol. Seventy-two hours after the cells were infected with lentivirus, 2 μg/mL puromycin was added to kill the cells that had not been transfected. shRNA sequences are provided (see Additional file 1: Table S2).

Cell counting kit-8 assay
U87 and U251 cells were assessed with the CCK-8 (Biosharp, China) reagent according to the manufacturer's instructions. Cells were inoculated on 96-well plates at a density of 2000 cells per well with 100 μL of medium. CCK8 solution (10 μL) was added to each well every 24 h for a total of 96 h, and the cells were further incubated at 37 °C for 1 h. The absorbance of each well was measured at 450 nm with a spectrophotometer.

Colony formation assay
U87 and U251 cells were prepared into a single cell suspension and seeded into a six-well plate (200 cells/well) for a two-week incubation to form colonies. After staining with 0.01% crystal violet (Sigma), the colonies were subjected to microscopic examination. The rate of colony formation was calculated.

Cell invasion and migration assays
After starving the cells for 6-8 h in serum-free DMEM, a total of 1 × 10 4 cells were seeded in the upper chamber with 200 μl of serum-free medium for the migration assay. In addition, 2 × 10 4 cells were added into Matrigel-coated upper transwell chambers for the invasion assay. The lower chambers were filled with DMEM containing 10% FBS. After incubation at 37 °C for 24 h, cells on the lower surface of the membrane were fixed in 100% methanol and stained with 0.1% crystal violet dye for 20 min at room temperature. Finally, after washing with phosphate-buffered saline, the cells were imaged in five randomly selected fields under a light microscope (Olympus Corporation) at × 100 magnification.

Flow cytometry cell cycle assay
After transient transfection, U87 and U251 cells were fixed in 75% ethanol for 12 h. Subsequently, cells were stained with propidium iodide (Beyotime) for cell cycle analysis. Finally, the percentage of cells in each cell cycle phase (G0/G1, S, and G2/M) was assessed, and the results were analyzed using the ModFit LT software.

RNA velocity and cells communication
The RNA velocity of the tumor cells was calculated using the package 'velocity' and 'scVelo' in Python. The various states of the GBM cells was mapped to show their internal transformation. The cross-talk between immunocytes and GBM cells was analyzed using the R package 'celltalker, ' and differential ligand-receptor pairs were identified.

Transcription factor (TF) regulatory network construction
The RcisTarget human database was downloaded from https:// resou rces. aerts lab. org/ cista rget/ for transcription factor regulatory network construction. The corresponding gene ranking motif database (Hg38_refseq-r80_10kb_ up_and_down_tss.mc9nr.feather, annotations_fname motifs-v9-nr.hgnc-m0.001-00.0.tbl) were downloaded from the human transcription factors list (https:/github. com/aertslab/pySCENIC/tree/master/resources), which is based on psSCENIC transcription factor regulation network. The AUCell algorithm was used to calculate the activity of each transcription factor, and the regulation module was identified according to the Connection Specificity Index (CSI). The calculation method of CSI was based on a previous article [17]. Similarly, we used the hTFtarget database to predict between TF and targets, which contains the most comprehensive data on human TF-target to date [18]. The overall activity score of the regulatory module was defined as the mean of all TF activities in the module.

Prediction of potential drug sensitivity
The drug sensitivity information and corresponding expression level were obtained from Genomics of Drug Sensitivity in Cancer (GDSC), Cancer Cell Line Encyclopedia (CCLE), and the Cancer Therapeutics Response Portal (CTRP) (https:// porta ls. broad insti tute. org/ ctrp). The CuAS score of each cell line was calculated and grouped based on the median. The correlation between the AUC and IC50 data of multiple drugs in the cell lines was calculated by using Spearman's correlation. The difference of the AUC value between the two groups were compared by the Wilcoxon test.

Statistical analysis
The significance of the difference between the two groups of continuous variables was evaluated using the Wilcox rank-sum test. Spearman's rank correlation was used to evaluate the correlation between the variables. Univariate and multivariate Cox regression and LASSO Cox regression were used to identify molecules with prognostic efficacy, and the K-M curves and log-rank tests were used to assess the survival differences between the sample groups. All computational analyses were performed by R (version 4.1.2) or Python.

Cuproptosis characteristic gene consistent clustering to identify sample subgroups
The RNA-seq data of 169 TCGA-GBM samples were obtained, and the tumor samples were clustered into two groups based on cuproptosis genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) through consistent cluster analysis (Fig. 1A, B, C). Significant differences were not observed in the survival between the two subgroups of samples, suggesting that these 10 genes alone may not be able to characterize the effect of cuproptosis mechanism on patient survival benefit (Fig. 1D). We also observed the expression patterns of these 10 characteristic genes in two subgroups of samples (Fig. 1E). However, significant differences were observed in the landscape of mutation, immune checkpoint expression level, and cancer hallmarks between the two subgroups. First, the SNV mutation frequencies of TP53 and other genes showed significant differences between the two types of samples ( Fig. 2A). The map of CNA showed that both types of patients had significant amplification on chromosome 7 (Fig. 2B), and significant difference was not observed in total frequency of CNA (Fig. 2C).
In addition, significant differences were observed in intratumoral heterogeneity between the two subgroups ( Fig. 2D). Moreover, we observed a significantly different expression level of immune checkpoint genes PD-1, IDO1, and LAG3 in the two subgroups (Fig. 2E). Significant differences were observed in the cancer hallmarks of fatty acid metabolism, KRAS, P53, NOTCH, and PI3K/ AKT/MTOR signaling pathway between the two subgroups (Fig. 2F).

Construction of cuproptosis activation scoring model based on differentially expressed genes between the sample subgroups
First, the DEGs between the two groups of samples were identified based on DESeq2 (Fig. 3A). The functions of these DEGs were enriched in the cell cycle related processes, protein kinase activity, P53 signaling pathway, and TGF-βsignaling pathway (Fig. 3B, C). In the TCGA-GBM sample set, we identified 14 candidate prognostic marker genes that were significantly associated with the OS of patients based on univariate Cox regression (Fig. 3E) and further filtered the redundant factors using LASSO Cox to obtain 11 prognostic marker genes (Fig. 3D). Based on the PCA of these 11 genes, their contribution to principal components 1 and 2 were used as coefficients (Fig. 3F) to construct CuAS.

Epiregulin (EREG) was an oncogenic gene that can influence immunity and cuproptosis
The EREG mRNA expression levels were high in tissues and multiple glioma cell lines ( Fig. 5 A, B). We used Western blot (Fig. 5C) and IHC staining ( Fig. 5E) to detect the EREG protein expression levels in tumors and normal tissues, and we found that the protein expression levels of EREG in tumors were higher than that in normal tissues. Additionally, we found that the protein expression levels of EREG in glioma cell lines were higher than that in normal astrocytic cell lines (Fig. 5G). Subsequently, we constructed knockdown stable cell lines of EREG and verified the knockdown effects on mRNA (Fig. 5F) and protein levels (Fig. 5D). Functional experiments showed that EREG knockdown (KD) can significantly inhibit the proliferation detected by Edu exepriments (Fig. 7E), invasion ( Fig. 6D), migration (Fig. 6E), and colony forming ability ( Fig. 6C) of tumor cells. Additionally, flow cytometry cell cycle assays suggested that EREG KD significantly inhibited cell cycle progression from the G0/G1 phase to the S phase (Fig. 7F). To explore the relationship between EREG and immune infiltration, we detected the expression level of PDL1 in EREG-KD group and found that PDL1 also decreased (Fig. 6A). To explore the relationship between EREG and cuproptosis, we performed the different treatment gradients of Cu-Elesclomol(ES) (1:1) on U251 cell lines and found that cell viability decreased with increasing time, and the effect of ES-Cu required a specific concentration range (5-50 nM) (Fig. 7A, B). Subsequently, we treated tumor cells with the same concentration (30 nM) with ES-Cu, and observed cell viability at 0, 12, 24, 36, 48, 60, 72, 84, and 96 h (Fig. 7D). It   chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9  chr10  chr11  chr12  chr13  chr14  chr15  chr16  chr17  chr18  c  was found that the proliferation rate of the treated cells decreased significantly compared with the cells that were not treated with ES-Cu. Then, the same treatment was performed on the shEREG and shNC groups and found that the proliferation rate of the two treated groups decreased significantly when compared with the shEREG and shNC groups that were not treated with ES-Cu; however, the reduction rate of the proliferation of the shEREG group was higher than that of the shEREG group with ES-Cu treatment, indicating that EREG can influence cell proliferation by affecting the process of cuproptosis (Fig. 7C, D). Therefore, we detected the protein expression level of FDX1 in the shEREG and shNC groups. The results showed that FDX1, the core regulatory protein in cuproptosis, was down-regulated in the shEREG group (Fig. 6B). Based on the above results, we believe that EREG is an oncogenic gene that can affect immunity by influencing the expression level of PDL1 and is closely related to the process of cuproptosis.

Single cell transcriptome analysis of CuAS patterns
Based on the downloaded single cell data (GSE173278, 29339 cells,10X Genomics platform), R Package Seurat was used to process the data. The expression profile was transformed by Log10, and 2,000 highly mutated genes were identified based on the VST method. Subsequently, principal component analysis and dimensionality reduction visualization were performed using UMAP. As a strong batch effect was observed, Harmony was used for batch correction (see Additional file 2: Fig. S1). Followup analysis was conducted based on the corrected data. The default parameters were used for clustering, and the meaning reference of analysis results based on known markers (SingleR was BP and HPCA) for cell type annotation (Fig. 8A). The cells in the GBM samples were divided into 7 categories, three malignant cell (OLIG1 + malignant, n = 11637; VEGFA + malignant, n = 6446; CENPF + malignant, n = 5363), microglia (n = 3219), fibroblasts (n = 1020), endothelial cells (n = 919), and oligodendrocytes (n = 735) (Fig. 8C). CuAS was calculated by using the previous model coefficients, but several cuproptosis characteristic genes were not detected in single cell data, and the expression level of many characteristic genes was undetectable (see Fig. 8E and Additional file 2: Fig. S2). A small number of cells (approximately 2800 cells) with a high CuAS score accounted for less than 10% of the whole cells and were distributed in multiple cell subpopulations. Most of the subpopulations contained less than 5% of cells with high CuAS, and 1132 cells with high CuAS were present in VEGFA + malignant cells (hypergeometric test, p value < 0.05). Trajectory inference of tumor cells are depicted (see S Additional file 2: Fig. S3). Moreover, the ancestor clone was determined based on CNV in combination with the idea of clone evolution, so as to determine the evolutionary relationship between cells more accurately. Based on CNV, we found that the OLIG1 + malignant cell may be the ancestor clone (Fig. 8D). Furthermore, in order to better explain the functional role of CuAS at the single cell level, the functions of VEGFA + malignant cells were observed and functional enrichment analysis was conducted based on specific up-regulated genes obtained by differential expression analysis. Pathways are mainly enriched in pathways related to hypoxia and oxidative stress (Fig. 8B).

Differential activation of transcription factors between high and low CuAS
Annotated files of human transcription factors were obtained from the RcisTarget database and the list of human transcription factors were downloaded. Transcription factor regulatory network was constructed using pySCENIC. Subsequently, the AUCell algorithm was used to calculate the activity of each transcription factor, and according to the CSI between the different transcription factors, four regulatory modules were identified (Fig. 9A). Module score was performed for each cell sample. We explored the association between cell type and module score, which revealed that the score of Module1 was significantly higher in VEGFA + malignant cells, while the score of Module2 was significantly higher in the endothelial cell subset. The score of Module3 was significantly higher in the microglia cell subset, while the score of Module4 was higher in the CENPF + malignant cells and partial OLIG1 + malignant cells, which reflected the differences of TF activated by different malignant cell subsets (Fig. 9B, C). Similarly, the results of the hTFtarget database showed that transcription factors such as FOSL2, JUND, NFIC and PBX3 were highly active in the VEGFA + malignant subgroup (see Additional file 2: Fig.  S4). In particular, we previously observed that cells with high CuAS scores were concentrated in VEGFA + malignant subsets, showing the potential association between CuAS scores and Module1.  samples, and immune-related pathways were differently enriched between the high and low CuAS groups, including T cells, NK Cell, B cell signal, chemokine signal, cytokine interaction, and other pathways (Fig. 10A). GO enrichment also showed that the activation differentiation and proliferation of T cells, NK cell proliferation, cytotoxic reaction, and other characteristics were highly enriched when CuAS scores were high (Fig. 10B). In addition, we identified differentially enriched signatures between the high and low CuAS groups based on GSEA enrichment analysis and also captured immune reaction processes such as leukocyte adhesion migration and T cell activation (Fig. 10C), indicating that the mechanism of cuproptosis is closely related to immune reaction process. Further, we calculated the immune and stromal components by using ESTIMATE, and it was observed that CuAS was significantly positively correlated with the stromal, immune, and ESTIMATE scores, while it was significantly negatively correlated with tumor purity (Fig. 10D-G), which suggested the association between cuproptosis and immunity. By using GSVA to calculate immune cell infiltration, we found that cuproptosis was significantly correlated with various types of immune cell infiltration, including activated DC and NK cells. (Fig. 12H). In addition, CIBERSORT and xCell methods were used to calculate various immune cell infiltrates, which revealed similar results (see Additional file 2: Fig.  S5).

Specific cell communication was different between high and low CuAS groups
As the high CuAS cells were significantly enriched in VEGFA + malignant cells, we mainly analyzed the difference in communication and function between the VEGFA + malignant and other cells. Extensive cell communication was observed in each cell subpopulation (Fig. 11A). Furthermore, by distinguishing between the incoming and outgoing signals, we found that fibroblasts are the dominant signaler of outgoing signaling, and VEGFA + malignant cells are the signal receivers (Fig. 11B). Furthermore, we identified two patterns of cell subpopulations in outgoing signaling, in which VEGFA + malignant cells belonged to Pattern 2 and corresponding pathways included VEGF, FGF, CDAM, CD22, ADGRE5, and other malignant progression related pathways (Fig. 11C). Meanwhile, we analyzed the dominant signaling pathway of each cell and found that VEGFA + malignant cells are not only involved in the signaling pathway of VEGF, but also in the CD99 signaling pathway, which was not proposed in the nonnegative Matrix Factorization (NMF) analysis (Fig. 11D). The CD99 signaling pathway plays an important role in tumor progression and transendothelial migration of cancer cells. VEGF and CD99 signaling pathways were further analyzed, and it was found that VEGFA + malignant cells are the dominant signalers of VEGF signals, and the cell subsets that were affected are mainly the endothelial cells and fibroblasts, both of which are important components of angiogenesis (Fig. 11 E, G). In addition, VEGFA + malignant cells are the dominant signaler, receiver, and influencer of CD99 signaling pathways, indicating that CD99 signaling pathways can occur as feedback loops (Fig. 11 F, H). Meanwhile, endothelial and fibroblast cells are also affected by CD99 signaling pathways, suggesting that VEGEA + malignant cells can influence transvascular endothelial migration (Fig. 12).
The immunofluorescence detection of tissue samples with high and low CuAS showed that VEGFA and CD99 were also highly expressed in tissues with high CuAS. The results were opposite in tissues with low CuAS (Fig. 13E, F, G), which provided a new idea for the intervention of cuproptosis-related tumor cells.

CuAS is associated with prognosis of immunotherapy
Based on all the immunotherapy data searched, we observed the ability of the CuAS score in predicting the prognosis and efficacy in the immunotherapy cohort. Phs001493 (Renal cell carcinoma, Anti-PD1 therapy) and PRJEB23709_ipiPD1 (Melanoma,anti-CTLA4 & AMP; Anti-pd1 dual antibody therapy) were significantly associated with worse prognosis (see Additional file 2: Fig.  S6B and D). For a patient's Progression Free Survival (PFS), we found that NCT02684006 (kidney cancer, anti-PDL1 treatment) was significantly associated with worse prognosis (see Additional file 2: Fig. S6E). Therefore, high CuAS patients may benefit from immunotherapy.

Potential targeted drugs for high CuAS glioblastoma cells
The expression data of cell lines were extracted from three databases: GDSC, CCLE, and CTRP. A lower AUC value represents a higher sensitivity to drugs. Using the AUC data provided by these databases, multiple drugs with a negative correlation between the AUC and signature were found in GDSC, such as methotrexate, BMS -708163, YM201636, FR -180204, and NU − 7441 (Fig. 12A). A variety of drugs with positive correlations were also found, such as cyclopamine (Fig. 12B). These drugs showed significant differences in the AUC between the groups of high and low signature (Fig. 12 C, D). No drugs with an IC50 significantly correlated with signature were found in CCLE, while drugs with a significant AUC (KU -0063794, cytochalasin B, GDC − 0941, cabozantinib, CI − 976, SJ -172550, SGX − 523, BRD − K71935468, temozolomide, AT7867, BRD-K66532283, palmostatin B, GDC-0879, ETP-46464, and NVP-BEZ235) negatively correlated with signature were found in CTRP (Fig. 12E,  F). These drugs were found in the AUC values of the high and low signature groups were significantly different in CTRP (Fig. 12 J, H). Therefore, high CuAS samples are likely to be sensitive to these compounds, and these compounds may be novel treatment options for GBM.

Experimental validation of model genes
The genes expression levels in the model were detected by qRT-PCR, and the results showed that they were highly expressed in 20 pairs of tumor and normal tissues (UNCX, SLC6A3, AGAP2-AS1, LINC00968, PTX3 and SBSPON), while ITPRID1, DCST2, ETV3L, and ENSG00000261327 were down-regulated (Fig. 13A).
According to the corresponding PCR results, we divided the tissue samples into high and low CuAS groups. IHC staining was performed on UNCX, SLC6A3, and PTX3, and it was found that the protein expression of the high CuAS group was higher than that of the low CuAS group (Fig. 13B, C, D).

Discussion
Copper is an essential cofactor in all organisms; however, it is toxic for cells when concentrations of copper exceed thresholds maintained by an evolutionarily conserved homeostasis mechanism [19,20]. In fact, it is not known how excessive copper can induce cell death. However, the Broad Institute has currently identified a new mechanism that is different from known cell death: cuproptosis [8]. Cuproptosis is a kind of cell death that is dependent on mitochondrial respiration. Copper directly binds to lipoylated components of the tricarboxylic acid cycle. Afterwards, aggregation of these copper-bound, lipoylated mitochondrial proteins and subsequent Fe-S cluster protein loss trigger proteotoxic stress and a distinct form of cell death [19][20][21][22]. Cuproptosis is involved in cell death, and the Broad Institute paper suggests that drugs that inhibit mitochondrial respiration may be a strategy against disease [19][20][21][22]. In addition, many mitochondrial proteins have a high degree of respiration function in various cancer cells [23]. Thus, copper ion metal carriers may be a new method for cancer treatment.
To the best of our knowledge, this study was the first paper to comprehensively analyze the association between copper-induced cell death and GBM by combining scRNA-seq and bulk RNA-seq data. First, we identified two sample subgroups based on the characteristic genes of cuproptosis. We found that immune checkpoint genes (PD-1, IDO1 and LAG3) and cancer hallmarks (fatty acid metabolism, KRAS, P53, NOTCH, and PI3K/ AKT/MTOR signaling pathway) showed significant differences between the two subgroups. Immune checkpoint is a kind of immunosuppressive molecule, which can regulate the intensity and breadth of the immune response, to avoid the damage and destruction of normal tissues. In the process of tumor occurrence and development, immune checkpoint has become one of the main reasons for immune tolerance. Subsequently, we constructed CuAS based on the differential genes of subgroups, which contained 11 genes, including 8 coding genes and 3 non-coding genes. EREG was the gene with the largest contribution coefficient to the principal component, so we focused on EREG. EREG is a 19-kDa peptide hormone that belongs to the Epidermal Growth Factor (EGF) family of peptide hormones [24]. Epiregulin binds to the EGF receptor (EGFR/ErbB1) and ErbB4 (HER4) and stimulates signaling of ErbB2 (HER2/Neu) and ErbB3 (HER3) through ligand-induced heterodimerization with a cognate receptor [24]. EREG possesses a range of functions in both normal physiologic states as well as in pathologic conditions. EREG contributes to inflammation, wound healing, tissue repair, and oocyte maturation by regulating angiogenesis and vascular remodeling and by stimulating cell proliferation [24]. Deregulated EREG activity appears to contribute to the progression of a number of different malignancies, including cancers of the bladder, stomach, colon, breast,  regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of transcription from RNA polymerase II promoter in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress regulation of DNA−templated transcription in response to stress response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion response to metal ion positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress positive regulation of transcription from RNA polymerase II promoter in response to stress regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway regulation of apoptotic signaling pathway cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to chemical stress cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels cellular response to decreased oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to oxygen levels cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia cellular response to hypoxia gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis gliogenesis apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes apoptotic mitochondrial changes negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway negative regulation of apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway intrinsic apoptotic signaling pathway glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation glial cell differentiation response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels response to decreased oxygen levels cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to inorganic substance cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion cellular response to metal ion response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to oxidative stress response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia response to hypoxia release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria release of cytochrome c from mitochondria autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion autophagy of mitochondrion mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly mitochondrion disassembly regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of oxidative stress−induced intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway regulation of intrinsic apoptotic signaling pathway PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response PERK−mediated unfolded protein response organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly organelle disassembly regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization regulation of mitochondrion organization response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to steroid hormone response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide response to hydrogen peroxide   lung, head and neck, and liver [2,7,24]. EREG is also associated with imaging omics as an important prognostic gene and MRI parameters revealed that hemodynamic abnormalities were associated with the expression level of the mTOR-EGFR pathway in patients with GBM [25]. Rab27b promotes the proliferation of adjacent cells and radio-resistance of highly malignant GBM cells through EREG-mediated paracrine signaling after irradiation [26]. Furthermore, EREG activates the extracellular signalingrelated kinase/MAPK pathway in GBM, suggesting that the inhibition of the EREG-EGFR interaction may be a strategy for EREG-overexpressing patients with GBM [2]. In our study, we detected EREG mRNA expression and protein levels in tissues and multiple glioma cell lines. IHC staining revealed that the EREG protein expression in tumors was higher than that in normal tissues; the result of WB also showed similar results. Knockdown of EREG can inhibit the proliferation, invasion, and migration of tumor cells. EGFR and PDL1 expression of protein were down-regulated after knockdowning of EREG. Moreover, we explored if EREG could influence the process of cuproptosis. Cell vitality assay demonstrated that only the coexistence of Cucl 2 and ES can influence the cell vitality and that other metals had no effect. The effect of ES-Cu required a specific concentration range (5 nM-50 nM). shEREG can revert the cell vitality that is influenced by cuproptosis. Therefore, we detected the protein expression of FDX1 in the shEREG and shNC groups. The results showed that FDX1, the core regulatory protein in cuproptosis, was down-regulated in the shEREG group.
Combined with the single cell transcriptome, the model of cuproptosis was analyzed, and the GBM sample cells were divided into seven types, including three types of malignant cells (OLIG1 + malignant, VEGFA + malignant, and CENPF + malignant). OLIG1 and other oligodendrocyte markers were highly expressed in OLIG1 + malignant cells, which may be oligodendrocyte progenitor glioma mother cells. VEGFA, CHI3L1, and other angiogenesis related markers were highly expressed in VEGFA + malignant cells, which may have a strong ability to induce local angiogenesis and may be associated with invasion/metastasis. CENPF + , TOP2A, UBE2C, and other markers are associated with the cell cycle and may be mesenchymal glioma blasts, which may be associated with tumor proliferation/invasion [27,28]. Others types observed were microglia, fibroblasts, endothelial cells, and oligodendrocytes. High CuAS was found in VEGFA + malignant cells. Based on CNV [29], OLIG1 + malignant cells were the ancestor clones. The function of VEGFA + malignant cells demonstrated that the pathways were mainly enriched in those related to hypoxia and stress, which is also consistent with the fact that cuproptosis is mitochondriondependent programmed cell death. Activated cells with high CuAS scores based on differences between high and low CuAS transcription factors were concentrated in the VEGFA + malignant cell subpopulation, reflecting the potential association between CuAS scores and Mod-ule1. The VEGF and CD99 signaling pathways were significantly enriched in high CuAS cells. VEGF specifically binds to Fltl and KDR/Flkl on the surface of endothelial cells, resulting in a variety of biological effects [30]. VEGF is closely associated with angiogenesis and development [31]. VEGF plays an important role in all stages of tumor formation, inducing the production of a large number of proteolytic enzymes, reducing the basement membrane of the host blood vessels, weakening the barrier effect, increasing the permeability of blood vessels, promoting a large amount of fibrinogen exudation, and forming a new matrix necessary for tumor adhesion and migration [30,31]. Angiogenesis is determined by the growth and metastasis of solid tumors. VEGF degrades extracellular matrix by inducing endothelial cells to express protease, resulting in metastasis, proliferation, and angiogenesis [32]. CD99 is abnormally expressed in many different types of tumors, and plays an important role in the diagnosis, development, metastasis, and prognosis, mainly affecting the invasion and metastasis of tumor cells [33]. Immunofluorescent detection of tissue samples with high and low CuAS showed that VEGFA and CD99 were also highly expressed in tissues with high CuAS, and the results were opposite in tissues with low CuAS, which provided a new idea for us to intervene in cuproptosisrelated tumor cells.
Immunotherapy is essential in tumor treatment. Despite the lack of specific immune cohort verification for glioma, several other tumor immune cohorts have shown the possibility of treatment for patients with high CuAS. Considering that EREG may affect the expression (See figure on next page.) Fig. 12 Drug sensitivity between different signature groups of GDSC and CTRP. A Negative correlation between signature score in GDSC and drug AUC (P < 0.05); B Positive correlation between signature score in GDSC and drug AUC (P < 0.05); C Differences in signature scores of GDSC cell lines with significant negative correlation under different drug treatments; D Differences in signature scores of cell lines with significant positive correlation under different drug treatments in GDSC. E Negative correlation between signature score and DRUG AUC in CTRP (P < 0.05); F Positive correlation between signature score in CTRP and drug AUC (P < 0.05); G Difference in signature scores of all cell lines in CTRP with significant negative correlation under different drug treatments; H Differences in signature scores of all cell lines in CTRP with significant positive correlation under different drug treatments